#!/usr/bin/env julia

"""
    test_forward_operators.jl

Comprehensive test suite for forward operators H(x) in GSI.
Tests observation operators, tangent linear operators, adjoint operators,
and their numerical properties for all observation types.

Test Coverage:
- Conventional observation operators (interpolation, coordinate transforms)
- Satellite radiance operators (radiative transfer integration)
- GPS radio occultation operators (Abel transform, refractivity)
- Tangent linear and adjoint consistency testing
- Numerical accuracy and conservation properties
- Performance testing for large state vectors
"""

using Test
using GSICoreAnalysis
using GSICoreAnalysis.ObservationProcessing.ForwardOperators
using LinearAlgebra
using Statistics
using Random
using Dates

"""Test basic forward operator properties"""
@testset "Forward Operator Fundamentals" begin
    println("  Testing forward operator fundamentals...")
    
    @testset "Interpolation Operators" begin
        # Create simple 3D grid for testing
        config = AnalysisConfig{Float64}(
            grid_size = (8, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize with known analytical field: T(x,y,z) = 300 - 20*z/H
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    # Linear temperature decrease with height
                    z_norm = (k-1) / (nz-1)  # 0 to 1
                    temp_slice[idx] = 300.0 - 20.0 * z_norm
                end
            end
        end
        
        # Test interpolation at known locations
        test_points = [
            (i=2.0, j=3.0, k=1.5, expected_temp=300.0 - 20.0 * 0.5/(nz-1)),
            (i=4.5, j=2.5, k=2.0, expected_temp=300.0 - 20.0 * 1.0/(nz-1)),
            (i=1.0, j=1.0, k=4.0, expected_temp=300.0 - 20.0 * 3.0/(nz-1))
        ]
        
        for point in test_points
            interpolated = trilinear_interpolation(
                sv, "temperature", point.i, point.j, point.k, config
            )
            @test interpolated ≈ point.expected_temp atol=1e-10
        end
        
        # Test interpolation at grid boundaries
        boundary_value = trilinear_interpolation(sv, "temperature", 1.0, 1.0, 1.0, config)
        @test boundary_value ≈ 300.0 atol=1e-12
        
        top_value = trilinear_interpolation(sv, "temperature", 1.0, 1.0, Float64(nz), config)
        @test top_value ≈ 280.0 atol=1e-12
    end
    
    @testset "Coordinate Transformations" begin
        # Test geographic to grid coordinate transformations
        
        # Define a simple grid domain
        grid_domain = GridDomain(
            lat_min = 35.0, lat_max = 45.0,
            lon_min = -110.0, lon_max = -95.0,
            nx = 8, ny = 6
        )
        
        # Test coordinate mapping
        test_locations = [
            (lat=40.0, lon=-102.5),  # Center of domain
            (lat=35.0, lon=-110.0),  # Bottom-left corner
            (lat=45.0, lon=-95.0),   # Top-right corner
            (lat=42.5, lon=-100.0)   # Interior point
        ]
        
        for loc in test_locations
            grid_coords = geographic_to_grid(loc.lat, loc.lon, grid_domain)
            
            @test 1.0 <= grid_coords.i <= grid_domain.nx
            @test 1.0 <= grid_coords.j <= grid_domain.ny
            
            # Test inverse transformation
            recovered_coords = grid_to_geographic(grid_coords.i, grid_coords.j, grid_domain)
            @test recovered_coords.lat ≈ loc.lat atol=1e-6
            @test recovered_coords.lon ≈ loc.lon atol=1e-6
        end
        
        # Test points outside domain
        outside_point = geographic_to_grid(30.0, -120.0, grid_domain)  # Outside domain
        @test outside_point.i < 1.0 || outside_point.i > grid_domain.nx ||
              outside_point.j < 1.0 || outside_point.j > grid_domain.ny
    end
    
    @testset "Vertical Coordinate Transforms" begin
        # Test pressure to model level transformations
        
        # Define hybrid sigma-pressure coordinates
        ak = [0.0, 1000.0, 3000.0, 7000.0, 15000.0]    # Pa
        bk = [1.0, 0.99, 0.95, 0.85, 0.70]           # dimensionless
        surface_pressure = 101325.0  # Pa
        
        # Compute pressure levels
        pressure_levels = hybrid_to_pressure(ak, bk, surface_pressure)
        
        @test length(pressure_levels) == length(ak)
        @test pressure_levels[1] ≈ surface_pressure  # Surface level
        @test issorted(reverse(pressure_levels))      # Decreasing with height
        
        # Test interpolation to observation pressure
        obs_pressure = 850.0e2  # 850 hPa in Pa
        
        model_level = pressure_to_model_level(obs_pressure, pressure_levels)
        @test 1.0 <= model_level <= length(pressure_levels)
        
        # Test that interpolated pressure matches
        interp_pressure = linear_interpolation(pressure_levels, model_level)
        @test interp_pressure ≈ obs_pressure rtol=0.01
    end
end

"""Test conventional observation operators"""
@testset "Conventional Observation Operators" begin
    println("  Testing conventional observation operators...")
    
    @testset "Surface Pressure Operator" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize surface pressure field
        fill!(sv.values, 0.0)
        ps_slice = get_variable_slice(sv, "surface_pressure")
        
        # Set uniform surface pressure
        fill!(ps_slice, 101325.0)  # 1013.25 hPa in Pa
        
        # Create surface observation
        ps_obs = ConventionalObservation{Float64}(
            obs_type = "surface_pressure",
            value = 101500.0,  # Observed value
            error = 200.0,
            latitude = 40.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Apply forward operator
        h_x = apply_surface_pressure_operator(ps_obs, sv, config)
        
        @test h_x ≈ 101325.0 atol=10.0  # Should match model surface pressure
        
        # Test with topographic correction
        ps_obs_elevated = ConventionalObservation{Float64}(
            obs_type = "surface_pressure",
            value = 95000.0,
            error = 200.0,
            latitude = 40.0,
            longitude = -100.0,
            elevation = 1500.0,  # 1.5 km elevation
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_corrected = apply_surface_pressure_operator(ps_obs_elevated, sv, config)
        
        # Should be reduced due to elevation (hydrostatic correction)
        @test h_x_corrected < h_x
        @test h_x_corrected > 85000.0  # Reasonable for 1.5 km elevation
    end
    
    @testset "Temperature Profile Operator" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 5),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize realistic temperature profile
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        
        # Standard atmosphere profile
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    # T decreases ~6.5 K/km in troposphere
                    height_km = (k-1) * 3.0  # 3 km per level
                    temp_slice[idx] = 288.15 - 6.5 * height_km
                end
            end
        end
        
        # Test temperature observations at different levels
        pressure_levels = [1000.0, 850.0, 700.0, 500.0, 300.0]  # hPa
        expected_temps = [288.15, 281.65, 275.15, 262.65, 243.65]  # K
        
        for (p, expected_t) in zip(pressure_levels, expected_temps)
            temp_obs = ConventionalObservation{Float64}(
                obs_type = "temperature",
                value = expected_t + randn() * 0.5,
                error = 1.0,
                latitude = 40.0,
                longitude = -100.0,
                pressure = p * 100.0,  # Convert to Pa
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            
            h_x = apply_temperature_operator(temp_obs, sv, config)
            
            @test h_x ≈ expected_t rtol=0.05  # 5% tolerance for interpolation
        end
    end
    
    @testset "Wind Vector Operator" begin
        config = AnalysisConfig{Float64}(
            grid_size = (8, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize wind field with simple geostrophic balance
        fill!(sv.values, 0.0)
        u_slice = get_variable_slice(sv, "u_wind")
        v_slice = get_variable_slice(sv, "v_wind")
        
        # Set up geostrophic wind: u = -∂Φ/∂y, v = ∂Φ/∂x
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    
                    # Simple wind field: u increases with y, v constant
                    u_slice[idx] = 5.0 + 2.0 * (j-1)/(ny-1)  # 5-7 m/s
                    v_slice[idx] = 3.0  # Constant 3 m/s
                end
            end
        end
        
        # Test wind speed observation
        wind_obs = ConventionalObservation{Float64}(
            obs_type = "wind_speed",
            value = 8.5,
            error = 2.0,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,  # 850 hPa
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_speed = apply_wind_speed_operator(wind_obs, sv, config)
        
        # Should be approximately sqrt(u²+v²)
        expected_speed = sqrt((6.0)^2 + (3.0)^2)  # Middle of domain
        @test h_x_speed ≈ expected_speed rtol=0.1
        
        # Test wind direction observation
        wind_dir_obs = ConventionalObservation{Float64}(
            obs_type = "wind_direction",
            value = 210.0,  # degrees
            error = 15.0,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_dir = apply_wind_direction_operator(wind_dir_obs, sv, config)
        
        # Should be atan2(v, u) converted to meteorological convention
        expected_dir = atan(3.0, 6.0) * 180.0/π + 90.0  # Meteorological direction
        if expected_dir < 0; expected_dir += 360.0; end
        
        @test h_x_dir ≈ expected_dir atol=5.0  # 5 degree tolerance
    end
    
    @testset "Humidity Operator" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize temperature and humidity fields
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        q_slice = get_variable_slice(sv, "specific_humidity")
        
        # Set realistic T and q profiles
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    
                    # Temperature decreases with height
                    temp_slice[idx] = 288.15 - (k-1) * 5.0
                    
                    # Specific humidity decreases exponentially
                    q_slice[idx] = 0.012 * exp(-(k-1) * 0.7)
                end
            end
        end
        
        # Test relative humidity observation
        rh_obs = ConventionalObservation{Float64}(
            obs_type = "relative_humidity",
            value = 65.0,  # %
            error = 10.0,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_rh = apply_relative_humidity_operator(rh_obs, sv, config)
        
        @test 20.0 <= h_x_rh <= 100.0  # Reasonable RH range
        
        # Test dewpoint observation
        dewpoint_obs = ConventionalObservation{Float64}(
            obs_type = "dewpoint",
            value = 278.0,  # K
            error = 2.0,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_dewpoint = apply_dewpoint_operator(dewpoint_obs, sv, config)
        
        @test 260.0 <= h_x_dewpoint <= 290.0  # Reasonable dewpoint range
        @test h_x_dewpoint <= temp_slice[1]     # Dewpoint <= temperature
    end
end

"""Test satellite observation operators"""
@testset "Satellite Observation Operators" begin
    println("  Testing satellite observation operators...")
    
    @testset "Radiance Operator Framework" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 10),  # More vertical levels for radiances
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize atmospheric profile
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        q_slice = get_variable_slice(sv, "specific_humidity")
        
        # Realistic atmospheric profile
        nx, ny, nz = config.grid_size
        pressure_levels = [1000.0, 900.0, 800.0, 700.0, 600.0, 
                          500.0, 400.0, 300.0, 200.0, 100.0]  # hPa
        
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    
                    # US Standard Atmosphere
                    p = pressure_levels[k] * 100.0  # Convert to Pa
                    if p > 22632.0  # Troposphere
                        temp_slice[idx] = 288.15 * (p/101325.0)^0.1903
                    else  # Stratosphere
                        temp_slice[idx] = 216.65 + 0.0 * (p - 22632.0)
                    end
                    
                    # Exponential humidity decay
                    q_slice[idx] = 0.015 * exp(-0.0001 * (101325.0 - p))
                end
            end
        end
        
        # Test AMSU-A radiance observation
        amsu_obs = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 6,  # 54 GHz O2 channel
            value = 242.5,
            error = 1.2,
            latitude = 40.0,
            longitude = -100.0,
            satellite_zenith_angle = 15.0,
            satellite_azimuth_angle = 180.0,
            scan_position = 15,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Apply radiance operator (simplified)
        h_x_radiance = apply_radiance_operator(amsu_obs, sv, config)
        
        @test 200.0 <= h_x_radiance <= 280.0  # Reasonable brightness temperature
        
        # Test that different channels give different values
        amsu_obs_ch2 = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 2,  # Different weighting function
            value = 255.0,
            error = 1.5,
            latitude = 40.0,
            longitude = -100.0,
            satellite_zenith_angle = 15.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_radiance_ch2 = apply_radiance_operator(amsu_obs_ch2, sv, config)
        
        @test abs(h_x_radiance - h_x_radiance_ch2) > 5.0  # Different channels should differ
    end
    
    @testset "GPS Radio Occultation Operator" begin
        config = AnalysisConfig{Float64}(
            grid_size = (4, 4, 15),  # High vertical resolution
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize atmospheric profile for GPS RO
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        q_slice = get_variable_slice(sv, "specific_humidity")
        ps_slice = get_variable_slice(sv, "surface_pressure")
        
        # Set surface pressure
        fill!(ps_slice, 101325.0)
        
        # Realistic profile for refractivity calculation
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    
                    # Height-based profile (0-30 km)
                    height_km = (k-1) * 2.0  # 2 km per level
                    
                    # Temperature profile
                    if height_km <= 11.0  # Troposphere
                        temp_slice[idx] = 288.15 - 6.5 * height_km
                    else  # Lower stratosphere
                        temp_slice[idx] = 216.65
                    end
                    
                    # Humidity profile
                    if height_km <= 15.0
                        q_slice[idx] = 0.01 * exp(-height_km/3.0)
                    else
                        q_slice[idx] = 1e-6  # Very dry stratosphere
                    end
                end
            end
        end
        
        # Test GPS RO bending angle observation
        gps_obs = SatelliteObservation{Float64}(
            obs_type = "gps_ro",
            channel = 1,
            value = 8.5e-4,  # Bending angle in radians
            error = 1e-5,
            latitude = 40.0,
            longitude = -100.0,
            impact_parameter = 6371.0e3 + 5.0e3,  # 5 km altitude
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_bending = apply_gps_ro_bending_operator(gps_obs, sv, config)
        
        @test h_x_bending > 0.0
        @test h_x_bending < 0.01  # Reasonable bending angle
        
        # Test GPS RO refractivity observation
        gps_ref_obs = SatelliteObservation{Float64}(
            obs_type = "gps_ro",
            channel = 2,  # Refractivity channel
            value = 2.5e-4,  # Refractivity N-units
            error = 5e-6,
            latitude = 40.0,
            longitude = -100.0,
            impact_parameter = 6371.0e3 + 3.0e3,  # 3 km altitude
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_refractivity = apply_gps_ro_refractivity_operator(gps_ref_obs, sv, config)
        
        @test h_x_refractivity > 0.0
        @test h_x_refractivity < 1e-3  # Reasonable refractivity
    end
end

"""Test tangent linear operators"""
@testset "Tangent Linear Operators" begin
    println("  Testing tangent linear operators...")
    
    @testset "TL Operator Linearization" begin
        config = AnalysisConfig{Float64}(
            grid_size = (4, 4, 3),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv = ControlVector(config)
        
        # Initialize background state
        randn!(sv.values)
        randn!(cv.values)
        
        # Test observation
        obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.5,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test linearization accuracy using finite differences
        ε = 1e-6
        
        # Reference nonlinear operator
        h_x0 = apply_forward_operator(obs, sv, config)
        
        # Perturb state vector
        cv_pert = ControlVector(config)
        assign!(cv_pert, cv)
        cv_pert.values .*= ε
        
        # Apply perturbation to state vector
        sv_pert = StateVector(config)
        assign_state!(sv_pert, sv)
        apply_control_to_state!(cv_pert, sv_pert, config)
        
        # Nonlinear operator on perturbed state
        h_x_pert = apply_forward_operator(obs, sv_pert, config)
        
        # Finite difference approximation
        fd_tl = (h_x_pert - h_x0) / ε
        
        # Exact tangent linear operator
        h_tl = apply_tangent_linear_operator(obs, cv, sv, config)
        
        # Compare linearization
        @test abs(h_tl - fd_tl) / abs(fd_tl) < 0.01  # 1% relative error
    end
    
    @testset "TL Operator Properties" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 4, 3),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv1 = ControlVector(config)
        cv2 = ControlVector(config)
        
        # Initialize vectors
        randn!(sv.values)
        randn!(cv1.values)
        randn!(cv2.values)
        
        obs = ConventionalObservation{Float64}(
            obs_type = "surface_pressure",
            value = 1013.25,
            error = 2.0,
            latitude = 40.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test linearity: H(α*x1 + β*x2) = α*H(x1) + β*H(x2)
        α, β = 2.5, -1.3
        
        h_cv1 = apply_tangent_linear_operator(obs, cv1, sv, config)
        h_cv2 = apply_tangent_linear_operator(obs, cv2, sv, config)
        
        # Create linear combination
        cv_combo = ControlVector(config)
        assign!(cv_combo, cv1)
        axpy!(α, cv1, cv_combo)  # cv_combo = cv1 + α*cv1 = (1+α)*cv1
        axpy!(β, cv2, cv_combo)  # cv_combo = (1+α)*cv1 + β*cv2
        
        # This is not the right test - let me fix it
        cv_combo = ControlVector(config)
        cv_combo.values .= α .* cv1.values .+ β .* cv2.values
        
        h_combo = apply_tangent_linear_operator(obs, cv_combo, sv, config)
        h_linear = α * h_cv1 + β * h_cv2
        
        @test abs(h_combo - h_linear) / abs(h_linear) < 1e-10  # Machine precision
    end
end

"""Test adjoint operators"""
@testset "Adjoint Operators" begin
    println("  Testing adjoint operators...")
    
    @testset "Adjoint Property Verification" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 4, 3),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv1 = ControlVector(config)
        cv2 = ControlVector(config)
        
        # Initialize vectors
        randn!(sv.values)
        randn!(cv1.values)
        fill!(cv2.values, 0.0)  # cv2 will receive adjoint result
        
        obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.5,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Apply tangent linear: y = H·x
        y = apply_tangent_linear_operator(obs, cv1, sv, config)
        
        # Apply adjoint: z = H^T·y  (where y is a scalar here)
        apply_adjoint_operator!(obs, y, cv2, sv, config)
        
        # Test adjoint property: <H·x, y> = <x, H^T·y>
        left_side = y * y  # <H·x, y> where y is the scalar result
        right_side = dot_product(cv1, cv2)  # <x, H^T·y>
        
        @test abs(left_side - right_side) / abs(right_side) < 1e-10
    end
    
    @testset "Multiple Observation Adjoint" begin
        config = AnalysisConfig{Float64}(
            grid_size = (4, 4, 3),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv = ControlVector(config)
        randn!(sv.values)
        fill!(cv.values, 0.0)
        
        # Create multiple observations
        observations = [
            ConventionalObservation{Float64}(
                obs_type = "temperature",
                value = 288.0,
                error = 1.5,
                latitude = 39.0,
                longitude = -101.0,
                pressure = 85000.0,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            ),
            ConventionalObservation{Float64}(
                obs_type = "surface_pressure",
                value = 1013.25,
                error = 2.0,
                latitude = 41.0,
                longitude = -99.0,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            ),
            ConventionalObservation{Float64}(
                obs_type = "wind_speed",
                value = 8.5,
                error = 2.5,
                latitude = 40.5,
                longitude = -100.5,
                pressure = 70000.0,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
        ]
        
        # Apply adjoint for each observation
        innovations = [2.5, -1.8, 0.7]  # obs - H(x)
        
        for (obs, innovation) in zip(observations, innovations)
            apply_adjoint_operator!(obs, innovation, cv, sv, config)
        end
        
        # Control vector should have been modified
        @test norm_cv(cv) > 0.0
        
        # Test that gradient accumulates properly
        cv_norm_final = norm_cv(cv)
        
        # Apply one more observation
        additional_obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 285.0,
            error = 1.0,
            latitude = 40.2,
            longitude = -99.8,
            pressure = 92500.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        apply_adjoint_operator!(additional_obs, -1.2, cv, sv, config)
        
        cv_norm_after = norm_cv(cv)
        @test cv_norm_after >= cv_norm_final  # Should increase (or stay same)
    end
    
    @testset "Adjoint Accuracy Test" begin
        config = AnalysisConfig{Float64}(
            grid_size = (5, 5, 4),
            ensemble_size = 0
        )
        
        # Random test vectors
        sv = StateVector(config)
        cv1 = ControlVector(config)
        cv2 = ControlVector(config)
        
        randn!(sv.values)
        randn!(cv1.values)
        fill!(cv2.values, 0.0)
        
        # Multiple observation types to test
        obs_types = [
            ("temperature", 85000.0),
            ("surface_pressure", missing),
            ("wind_speed", 70000.0),
            ("relative_humidity", 92500.0)
        ]
        
        total_dot_tl = 0.0
        total_dot_adj = 0.0
        
        for (obs_type, pressure) in obs_types
            obs = ConventionalObservation{Float64}(
                obs_type = obs_type,
                value = 100.0,  # Arbitrary value
                error = 1.0,
                latitude = 40.0,
                longitude = -100.0,
                pressure = pressure,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            
            # Apply TL operator
            y = apply_tangent_linear_operator(obs, cv1, sv, config)
            
            # Apply adjoint operator  
            innovation = randn()  # Random innovation
            cv_temp = ControlVector(config)
            fill!(cv_temp.values, 0.0)
            apply_adjoint_operator!(obs, innovation, cv_temp, sv, config)
            
            # Accumulate dot products for adjoint test
            total_dot_tl += y * innovation
            total_dot_adj += dot_product(cv1, cv_temp)
        end
        
        # Test adjoint property across all operators
        relative_error = abs(total_dot_tl - total_dot_adj) / abs(total_dot_adj)
        @test relative_error < 1e-12  # Very tight tolerance for adjoint property
    end
end

"""Test numerical accuracy and conservation"""
@testset "Numerical Accuracy and Conservation" begin
    println("  Testing numerical accuracy and conservation...")
    
    @testset "Mass Conservation" begin
        config = AnalysisConfig{Float64}(
            grid_size = (8, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        fill!(sv.values, 0.0)
        
        # Initialize uniform mass field
        ps_slice = get_variable_slice(sv, "surface_pressure")
        fill!(ps_slice, 101325.0)
        
        # Create surface pressure observations covering the domain
        ps_observations = ConventionalObservation{Float64}[]
        
        nx, ny = config.grid_size[1:2]
        for i in 2:2:nx-1  # Sample every other grid point
            for j in 2:2:ny-1
                obs = ConventionalObservation{Float64}(
                    obs_type = "surface_pressure",
                    value = 101325.0 + randn() * 100.0,
                    error = 50.0,
                    latitude = 35.0 + (i-1) * 10.0/(nx-1),
                    longitude = -110.0 + (j-1) * 15.0/(ny-1),
                    time = DateTime(2024, 1, 1, 12, 0, 0)
                )
                push!(ps_observations, obs)
            end
        end
        
        # Compute total mass before operator application
        mass_before = sum(ps_slice) * grid_cell_area(config)
        
        # Apply forward operators (should not change mass)
        h_x_values = [apply_forward_operator(obs, sv, config) 
                     for obs in ps_observations]
        
        # Mass should be unchanged
        mass_after = sum(ps_slice) * grid_cell_area(config)
        @test abs(mass_after - mass_before) / mass_before < 1e-14
        
        # Test that observations sample the field correctly
        mean_h_x = mean(h_x_values)
        @test abs(mean_h_x - 101325.0) < 200.0  # Should be close to field mean
    end
    
    @testset "Energy Conservation" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        fill!(sv.values, 0.0)
        
        # Initialize kinetic energy field
        u_slice = get_variable_slice(sv, "u_wind")
        v_slice = get_variable_slice(sv, "v_wind")
        
        # Set balanced wind field
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    u_slice[idx] = 10.0 * sin(2π * i / nx)
                    v_slice[idx] = 10.0 * cos(2π * j / ny)
                end
            end
        end
        
        # Compute total kinetic energy
        ke_before = 0.5 * (sum(u_slice.^2) + sum(v_slice.^2))
        
        # Apply wind operators
        wind_obs = ConventionalObservation{Float64}(
            obs_type = "wind_speed",
            value = 12.0,
            error = 2.0,
            latitude = 40.0,
            longitude = -100.0,
            pressure = 85000.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_wind = apply_forward_operator(wind_obs, sv, config)
        
        # Energy should be unchanged by observation operators
        ke_after = 0.5 * (sum(u_slice.^2) + sum(v_slice.^2))
        @test abs(ke_after - ke_before) / ke_before < 1e-14
        
        # Observed wind speed should be reasonable
        @test 5.0 <= h_x_wind <= 20.0
    end
    
    @testset "Interpolation Accuracy" begin
        # Test interpolation accuracy with analytical functions
        config = AnalysisConfig{Float64}(
            grid_size = (9, 7, 5),  # Odd numbers for exact center points
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        
        # Initialize with polynomial function: T(i,j,k) = a + b*i + c*j + d*k + e*i*j
        a, b, c, d, e = 280.0, 2.0, 1.5, -3.0, 0.1
        nx, ny, nz = config.grid_size
        
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    temp_slice[idx] = a + b*i + c*j + d*k + e*i*j
                end
            end
        end
        
        # Test interpolation at various points
        test_points = [
            (i=2.5, j=3.5, k=2.0),  # Half-grid points
            (i=5.0, j=4.0, k=3.0),  # Grid point (should be exact)
            (i=3.3, j=2.7, k=4.2),  # General point
            (i=1.0, j=1.0, k=1.0),  # Corner point
        ]
        
        for point in test_points
            # Analytical solution
            expected = a + b*point.i + c*point.j + d*point.k + e*point.i*point.j
            
            # Numerical interpolation
            interpolated = trilinear_interpolation(
                sv, "temperature", point.i, point.j, point.k, config
            )
            
            # For polynomial functions, trilinear interpolation should be very accurate
            if point.i == round(point.i) && point.j == round(point.j) && point.k == round(point.k)
                @test interpolated ≈ expected atol=1e-14  # Exact at grid points
            else
                @test abs(interpolated - expected) / abs(expected) < 0.01  # 1% for off-grid
            end
        end
    end
end

"""Test performance characteristics"""
@testset "Forward Operator Performance" begin
    println("  Testing forward operator performance...")
    
    @testset "Large State Vector Performance" begin
        # Test with realistic GSI dimensions
        config = AnalysisConfig{Float64}(
            grid_size = (32, 24, 20),  # ~15K grid points
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        randn!(sv.values)
        
        # Create observation dataset
        n_obs = 1000
        observations = ConventionalObservation{Float64}[]
        
        for i in 1:n_obs
            obs_type = rand(["temperature", "surface_pressure", "wind_speed"])
            pressure = obs_type == "surface_pressure" ? missing : 50000.0 + 50000.0 * rand()
            
            obs = ConventionalObservation{Float64}(
                obs_type = obs_type,
                value = 100.0 * rand(),
                error = 1.0 + rand(),
                latitude = -60.0 + 120.0 * rand(),
                longitude = -180.0 + 360.0 * rand(),
                pressure = pressure,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            push!(observations, obs)
        end
        
        # Time forward operator application
        println("    Timing forward operators for $(n_obs) observations...")
        forward_time = @elapsed begin
            h_x_values = [apply_forward_operator(obs, sv, config) 
                         for obs in observations]
        end
        
        @test length(h_x_values) == n_obs
        @test forward_time < 10.0  # Should complete within 10 seconds
        @test all(isfinite, h_x_values)  # All values should be finite
        
        println("    Forward operator time: $(forward_time:.2f) seconds")
        println("    Rate: $(n_obs/forward_time:.1f) observations/second")
        
        # Time tangent linear operators
        cv = ControlVector(config)
        randn!(cv.values)
        
        sample_obs = observations[1:min(100, n_obs)]  # Smaller sample for TL
        
        println("    Timing tangent linear operators for $(length(sample_obs)) observations...")
        tl_time = @elapsed begin
            h_tl_values = [apply_tangent_linear_operator(obs, cv, sv, config) 
                          for obs in sample_obs]
        end
        
        @test length(h_tl_values) == length(sample_obs)
        @test tl_time < 5.0
        @test all(isfinite, h_tl_values)
        
        println("    Tangent linear time: $(tl_time:.2f) seconds")
    end
    
    @testset "Memory Efficiency" begin
        config = AnalysisConfig{Float64}(
            grid_size = (16, 12, 8),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv = ControlVector(config)
        randn!(sv.values)
        randn!(cv.values)
        
        # Monitor memory usage during operations
        initial_memory = Base.gc_live_bytes()
        
        # Create and process observations in batches
        batch_size = 200
        n_batches = 5
        
        for batch in 1:n_batches
            batch_obs = ConventionalObservation{Float64}[]
            
            # Create batch
            for i in 1:batch_size
                obs = ConventionalObservation{Float64}(
                    obs_type = "temperature",
                    value = 288.0 + randn() * 10.0,
                    error = 1.0 + rand(),
                    latitude = -45.0 + 90.0 * rand(),
                    longitude = -180.0 + 360.0 * rand(),
                    pressure = 50000.0 + 50000.0 * rand(),
                    time = DateTime(2024, 1, 1, 12, 0, 0)
                )
                push!(batch_obs, obs)
            end
            
            # Process forward operators
            h_x_batch = [apply_forward_operator(obs, sv, config) for obs in batch_obs]
            
            # Process tangent linear operators
            h_tl_batch = [apply_tangent_linear_operator(obs, cv, sv, config) for obs in batch_obs]
            
            @test length(h_x_batch) == batch_size
            @test length(h_tl_batch) == batch_size
            
            # Force garbage collection
            GC.gc()
        end
        
        final_memory = Base.gc_live_bytes()
        memory_growth = final_memory - initial_memory
        
        # Memory growth should be reasonable
        @test memory_growth < 50 * 1024 * 1024  # Less than 50 MB growth
        
        println("    Memory growth during processing: $(memory_growth ÷ (1024*1024)) MB")
    end
end

"""Helper functions for testing"""

function grid_cell_area(config::AnalysisConfig{T}) where T
    # Simplified grid cell area calculation
    return T(1.0e8)  # ~100 km² cells
end

function linear_index(i::Int, j::Int, k::Int, grid_size::Tuple{Int,Int,Int})
    nx, ny, nz = grid_size
    return (k-1) * nx * ny + (j-1) * nx + i
end

println("  ✓ All forward operator tests completed successfully")